#include "matrixMath.h"

namespace moremath{

	Matd getInnerProdMat(vector<Vecd> _a, vector<Vecd> _b){
		//assert(_a.size()==_b.size());
		Matd matG(_a.size(), _b.size(), vl_0);
		for(int i=0; i<_a.size(); i++){
			for(int j=0; j<_b.size(); j++){
				matG[i][j] = dot(_a[i], _b[j]);
			}
		}
		return matG;
	}

	Matd getInnerProdMat(Matd _a, Matd _b){
		Matd matG(_a.Cols(), _b.Cols(), vl_0);
		for(int i=0; i<matG.Rows(); i++){
			for(int j=0; j<matG.Cols(); j++){
				matG[i][j] = dot(col(_a, i), col(_b, j));
			}
		}
		return matG;
	}

	void getColMeanCenteredMat(Matd &_m, Matd *_newm, Vecd *_mean){
		*_newm = _m;
		// compute the mean;
		*_mean = Vecd(_m.Rows(), vl_0);
		for(int j=0; j<_m.Cols(); j++){
			*_mean += col(_m, j);
		}
		*_mean = *_mean/_m.Cols();
		// subtract the mean
		for(int j=0; j<_m.Cols(); j++){
			for(int i=0; i<_m.Rows(); i++){
				(*_newm)[i][j] = _m[i][j] - (*_mean)[i];
			}
		}
	}
	// trace of A*B
	double traceMult(Matd &_A, Matd &_B){
		double t = 0;
		for(int i=0; i<_A.Rows(); i++){
			for(int j=0; j<_A.Cols(); j++){
				t += _A[i][j]*_B[j][i];
			}
		}
		return t;
	}
	// trace of A*B^T
	double traceMultT(Matd &_A, Matd &_B){
		double t = 0;
		for(int i=0; i<_A.Rows(); i++){
			for(int j=0; j<_A.Cols(); j++){
				t += _A[i][j]*_B[i][j];	// j,ith element of B^T
			}
		}
		return t;
	}

}
